/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package exercicio.util;

/**
 * @author Herminio S. Gomes
 */
public class Distri {

    final static double r2 = 0.70710678118654752440;   // r2=raiz(2)/2
    
    /**
     * Calcula C(n, k).
     * @param n Possibilidades.
     * @param k Saídas.
     * @return A combinação C(n, k).
     */
    public static double comb(int n, int k) {
        // { boolean erro = (n<0) || (k<0) || (n<k); if (erro) System.exit(1);}
        double c = 1;
        int ka = (n - k < k) ? n - k : k;
        if (ka == 0) {
            return c;
        }
        for (int i = 1; i <= ka; i++) {
            c = c * (n - i + 1) / i;
        }
        return c;
    }
    
    /**
     * Calcula o fatorial de um número.
     * @param n Número.
     * @return Fatorial do número (n!).
     */
    public static double fatorial(int n) {
        if (n <= 1) {
            return 1;
        }
        double f = 1;
        for (int i = 2; i <= n; i++) {
            f = f * i;
        }
        return f;
    }

    public static double binomCDF(int n, double p, int k1, int k2) {
        //  { boolean erro = (n<0) || (k1<0) || (k2<0) || (n<k1)
        //    || (n<k2) || ( k1>k2); if (erro) System.exit(1);}
        double s = 0, q = 1 - p;
        for (int k = k1; k <= k2; k++) {
            s += comb(n, k) * Math.pow(p, k) * Math.pow(q, n - k);
        }
        return s;
    }

    public static double binomPDF(int n, double p, int k) {
        return binomCDF(n, p, k, k);
    }

    public static double hipergeomPDF(int k, int m, int N, int n) {
        // P(x=k)=comb(m,k)*comb(N-m,n-k)/comb(N,n)
        //  { boolean erro = (n<0) || (k<0) || (k>n) || (k>m)
        //       || (m>N) || (n>N) || ( N<2); if (erro) System.exit(1);}

        return comb(m, k) * comb(N - m, n - k) / comb(N, n);
    }

    public static double hipergeomCDF(int k1, int k2, int m, int N, int n) {
        // { boolean erro = (k2<k1); if(erro) System.exit(1);
        double s = 0;
        for (int k = k1; k <= k2; k++) {
            s += hipergeomPDF(k, m, N, n);
        }
        return s;
    }

    public static double funcaoErro(double x) {
        double ef, tt, t, xx,
                p = 0.3275911, a1 = 0.254829592, a2 = -0.284496736,
                a3 = 1.421413741, a4 = -1.453152027, a5 = 1.061405429;
        if (x > 0) {
            xx = x;
        } else {
            xx = -x;
        }
        t = 1 / (1 + p * xx);
        tt = ((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t;
        ef = 1 - Math.exp(-xx * xx) * tt;
        if (x < 0) {
            return -ef;
        }
        return ef;
    }

    public static double invErro(double r) {
        // { boolean erro = (Math.abs(r)>1); if (erro) System.exit(1);}
        return r2 * invnormal((r + 1) / 2.0);
    }

    public static double normalPDF(double z) {
        return Math.exp(-z * z / 2) / Math.sqrt(2 * Math.PI);
    }

    public static double normalCDF(double z) {
        if (z < 0) {
            return 1 - normalCDF(-z);
        }
        double b1 = 0.319381530, b2 = -0.356563782, b3 = 1.781477937,
                b4 = -1.821255978, b5 = 1.330274429, t = 1 / (1 + 0.2316419 * z);
        return 1 - normalPDF(z) * (((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t);
    }

    public static double normalCDF(double z1, double z2) {
        // { boolean erro = (z1>z2); if (erro) System.exit(1);}
        return normalCDF(z2) - normalCDF(z1);
    }

    public static double normalCDF(double x1, double x2, double media, double desvio) {
        // { boolean erro = (x1>x2) || (desvio<=0) ; if (erro) System.exit(1);}
        double z1 = (x1 - media) / desvio, z2 = (x2 - media) / desvio;
        return normalCDF(z1, z2);
    }

    public static double normalCDF(double x, double media, double desvio) {
        // { boolean erro = (desvio<=0); if (erro) System.exit(1);}
        return normalCDF((x - media) / desvio);
    }

    public static double invnormal(double p) {
        // { boolean erro = (p<0) || (p>1); if (erro) System.exit(1);}
        if (p > 0.999968) {
            return 4;
        }
        if (p < 0.000032) {
            return -4;
        }
        if (p > 0.5) {
            return -invnormal(1 - p);
        }
        double t = Math.sqrt(Math.log(1 / (p * p)));
        double c0 = 2.515517, c1 = 0.802853, c2 = 0.010328,
                d1 = 1.432788, d2 = 0.189269, d3 = 0.001308;
        double pn = c0 + t * (c1 + c2 * t), pd = ((d3 * t + d2) * t + d1) * t + 1;
        double x = pn / pd - t;
        // refina  x  fazendo tres iterações de Newton
        for (int j = 1; j <= 3; j++) {
            x = x - (normalCDF(x) - p) / normalPDF(x);
        }
        return x;
    }

    public static double exponencialCDF(double x, double media) {
        return 1 - Math.exp(-x / media);
    }

    public static double exponencialCDF(double x1, double x2, double media) {
        // { boolean erro = (x2<x1) || (media<=0); if (erro) System.exit(1);}
        return exponencialCDF(x2, media) - exponencialCDF(x1, media);
    }

    public static double poissonPDF(double lambda, int k) {
        return poissonCDF(lambda, k, k);
    }

    public static double poissonCDF(double lambda, int k1, int k2) {
        double p = 0, el = Math.exp(-lambda);
        for (int k = k1; k <= k2; k++) {
            p += Math.pow(lambda, k) / fatorial(k);
        }
        return p * el;
    }

    public static double rayleighPDF(double x, double sig) {
        double sig2 = sig * sig;
        return x / sig2 * Math.exp(-x * x / (2 * sig2));
    }

    public static double rayleighCDF(double x, double sig) {
        double sig2 = 2 * sig * sig;
        return 1 - Math.exp(-x * x / sig2);
    }

    public static double gamma(double x) {   // { boolean erro = (Math.abs(x) < 1e-4); if (erro) System.exit(1);}
        double p[] = {1.000000000190015, 76.18009172947146,
            -86.50532032941677, 24.01409824083091,
            -1.231739572450155, 1.208650973866179E-3, -5.395239384953E-6};
        double s = p[0], r2pi = 2.5066282746310005;
        for (int n = 1; n <= 6; n++) {
            s += p[n] / (x + n);
        }
        s *= r2pi / x;
        return s * Math.pow(x + 5.5, x + 0.5) * Math.exp(-x - 5.5);
    }

    public static double betaPDF(double a, double b, double x) {  //  Valor esperado =a /(a+b)
        // variancia = ab/( (a+b+1)(a+b)^2 )
        return beta(a, b) * Math.pow(x, a - 1) * Math.pow(1 - x, b - 1);
    }

    public static double beta(double a, double b) {
        return gamma(a + b) / (gamma(a) * gamma(b));
    }
}